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Abstract 

Kinetics of conformational change of a semiflexible polymer under mechanical external field were 
investigated with Langevin dynamics simulations. It is found that a semiflexible polymer exhibits 
large hysteresis in mechanical folding/unfolding cycle even with a slow operation, whereas in a 
flexible polymer, the hysteresis almost disappears at a sufficiently slow operation. This suggests 
that the essential features of the structural transition of a semiflexible polymer should be interpreted 
at least on a two-dimensional phase space. The appearance of such large hysteresis is discussed in 
relation to different pathways in the loading and unloading processes. By using a minimal two- 
variable model, the hysteresis loop is described in terms of different pathways on the transition 
between two stable states. 

PACS numbers: 36.20.-r,87.15.La,05.70.Fh 
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I. INTRODUCTION 



Recent developments in experimental techniques regarding single molecule observation 
and manipulation have enabled us to directly investigate the statistical properties of indi- 
vidual small systems. 

ELQ 

These experiments have provided novel insights into the structure 
and function of biological macromolecules. These studies have revealed significant effects of 
fluctuation on the properties of single molecules. 

Most biopolymers, including DNA and many proteins, have finite stiffness and are there- 
fore classified as semiflexible polymers. It is expected that the stiffness plays an important 
role in the structural stability and function of these biomacromolecules. Recent experiments 
and simulations have clarified that a single semiflexible polymer exhibits first order phase 
transition between elongated coil and folded compact states with a decrease in the solvent 
quality. A single semiflexible polymer folds into various kinds of ordered structures de- 
pending on its stiffness and temperature, such as a toroid or a rod. ^ This implies that the 
density of monomers is not sufficient to characterize ordered structures, but the orienta- 
tional order plays an essential role in the folding transition of a semiflexible chain, i.e., the 
density-order coupling should be crucial. 

Among the studies on single molecule manipulation, experiments on mechanical unfolding 
with laser tweezers or atomic force microscopy have been actively performed over the past 
decade. In these experiments, a non-trivial force response, which is called a sawtooth- 

like pattern or a stick- and-release pattern, has been observed either in a protein or in a 
DNA molecule. Some authors claim that such patterns are associated with repeated domain 
structures in a m o,ecule.aa In — to this, it is found in flQ that a DNA ntoiecule 
having no repeated sequence exhibits a stick- and-release pattern at the collapsing conditions 
where the concentration of multivalent cations is high, whereas the mechanical response of 
a worm-like chain is observed when a DNA molecule is in the coil conditions at the low 
concentration of multivalent cations. Furthermore, it has been argued theoretically that a 
toroidal structure of a folded semiflexible polymer without any repeated domains leads to a 



stick- and-release pattern. 
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A system under loading and unloading processes is clearly in a nonequilibrium state 
unless the rate of change is infinitesimally slow, and therefore the kinetics play essential 
roles. Recently, it has been revealed that a difference in the equilibrium free energy can 
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and the theory 



be calculated from kinetic processes under a nonequilibrium state 

has been applied to biopolymer stretching. Q] Despite such an analysis of nonequilibrium 
effects, there have been few studies on the large hysteresis loop observed in the loading and 
unloading cycle of biomolecules. 

Hysteresis is a measure of how the system is away from thermal equilibrium. At the 
length scale of micron or submicron, thermalization is so fast that equilibrium states are 
easily attained. However, to perform work against an external force, as seen with molecular 
motors in biological systems, the system should remain in a nonequilibrium state. Therefore, 
it is important to clarify the origin of hysteresis and, especially, its robustness. 

Theoretical and numerical studies are needed to gain general insight into not only the 
nonequilibrium characteristics of a single polymer chain, but also the working mechanism of 
a macromolecule in biological systems. In the present article, we consider flexible and semi- 
flexible polymers under strain by molecular dynamics simulations. We investigate the force 
response and focus on the close connection between hysteresis and the structural transition 
of a single polymer chain. 

A coarse-grained approach is generally quite useful to understand the essence of the 
phase transitions. We shall introduce a Ginzburg-Landau type free energy and the kinetic 
equations for the present problem. As mentioned above, however, the density of monomers 
is not sufficient to distinguish the different conformation in a semiflexible chain. Therefore, 
we introduce another kinetic variable corresponding to the orientational order. We do not 
intend to derive the free energy starting from the chain model but show that the coarse- 
grained free energy having two metastable states indeed describes properly the hysteresis in 
the force-response relation. To our knowledge, this kind of approach has not been available 
in the kinetics of a single chain. 

In the next section, we present our chain model and describe the method of simulations. 
In section III, the results of simulations are given. A phenomenological scaling argument 
is developed in section IV. The two-variable model and its numerical studies are shown in 
section V. A summary is given in section VI. 
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II. SIMULATIONS 
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We adopt the parameters to be 



To examine the folding and unfolding kinetics, we carried out Langevin dynamics simu- 
lations of a bead-spring model using the following potentials 

k a \ 

Vbeads = y 2j(|r i+ i - r»| - af (1) 

i 

Kend = |^(l-COS^) 2 (2) 

i 

VLj = 4e£«^S) 12 - (^^) 6 ), (3) 

^ Tj-rJ \Ti - rA 

where V = Vb ea ds + K>end + VLj and rj is the coordinate of the ith monomer and 6i is the 
angle between adjacent bond vectors. The monomer size a and ksT are chosen as the unit 
length and energy, respectively, where ks is the Boltzmann constant and T is the absolute 
temperature. Monomer-monomer interaction is included by the Lennard- Jones potential 
controlled by e. A similar bead-spring model JTJ) - (J3J) has been used to study the properties 
of semiflexible chains by setting a large value of k a . 
k a = 400 and e = 1.0. The bending elasticity is chosen to be k = 100 for a semiflexible 
polymer and k = for a flexible polymer. The persistence length l p is a convenient measure 
to characterize the stiffness of a polymer chain. In the present model, l p ~ 13.5a for a 
semiflexible polymer. Because of discritization of a chain with beads, the minimum of 
the persistence length is a, and therefore l p ~ a for a flexible polymer. We consider a 
homopolymer with a polymerization index N = 200. 

A remark is now in order. Because a semiflexible chain is locally inextensible, the partition 
function should be written as Z = j T>Y,jU.5(\Y i+ i — rj| — a)e _ ^ bcnd+VLJ ^. As an analytically 
or numerically tractable model, several approximations have been employed and examined. 
For example, the local constraint was replaced by a global constraint in refs. jl^l. 17] to 
study analytically the equilibrium properties of a semiflexible polymer. Instead of such 
approximations, we adapt the harmonic potential (JTJ) between a pair of the neighboring 
beads. A continuum version of the model Q was introduced to evaluate the end-to-end 
distance which shows clearly that the end-to-end distance is insensitive to k a for large 
values of k a . 

The equation of motion can be written as, 

d 2 ri dr; dV 



where m and 7 are the mass and friction constant of a monomeric unit, respectively. The 
constant r = 7a 2 //ceT is chosen to be the unit for the time scale. We set the time step 
as O.Olr, and use m = 1.0 and 7 = 1.0. Gaussian white noise £j satisfies a fluctuation- 
dissipation relation, 



where < ■ • • > indicates an ensemble average. 

The folding and unfolding states of a polymer are characterized by the average local 
monomer density p, which is defined by, 



where H(x) is the Heaviside function. In this work, we set r c = 3.0. 

Several comments on the values of the parameters chosen above are in order. The at- 
tractive strength e in the Lennard- Jones potential © is chosen so that ordered structures 
such as toroids and rods are formed for polymerization index N = 200 and l v ~ 13.5a. |^ 
Among such ordered structures, a toroidal structure is the most energetically favorable in 
the present set of parameters. The spring constant k a in is set to be a sufficiently large 
value to realize the local inextensibility of a chain. We set m = 7 = 1.0 in Q to opti- 
mize the numerical simulations within accessible simulation time avoiding any numerical 
instabilities. With this condition, the inertia term does not cause any artifacts because the 
relaxation time of the momentum of a monomer is sufficiently faster compared to that of 
conformational change. 

We fix one end of a polymer chain and move the other end at constant velocity z. The 
end-to-end distance z is one of the measurable quantities. The force / applied to the 
iVth monomer is monitored during the operation. Note that force is averaged over time 
to avoid its large fluctuation. The inverse velocity A = z~ x is often used as a measure of 
the speed of external operations. It is convenient to define Afoid, which corresponds to the 
time scale of the folding transition of a semiflexible polymer in bulk. It can be written as 
Afoid =< r foid/(-R C oii — -Rfoid) >) where Tfoid is the characteristic time for a folding transition, 
and i? CO ii and -Rfoid are the end-to-end distance of the coiled and folded states in a semiflexible 
polymer, respectively. 

Here we mention the time and length scales in our simulations. The advantage of the 
bead-spring model (jlj) with (JT]|3J) is that we can simulate on a realistic time scale. The folding 
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transition time is on the order of 10 5 — 10 6 steps, which corresponds to about 0.1-1.0 sec in 
the experiment with DNA. [3] From this correspondence, our loading and unloading speed 
(10 7 — 10 s steps in the simulations) is considered to correspond to 10 - 100 sec, which is a 
experimentally accessible value. We can also discuss the length scale. The persistence length 
of DNA is about 50 nm, which corresponds to about 10a in our simulation. Therefore, we can 
estimate the size of a toroidal state of a semiflexible polymer (~ 10a in the simulation) to be 
about 50 nm, which is in good agreement with experiments on DNA. Most simulations suffer 
from a large value of the force response. (^(J However, the force response in our simulations 
(2a — 5a/kBT) corresponds to 3 — 7pN, which is consistent with the experimental value. [3] 

III. RESULTS 

Figure Q shows the local monomer density and the force response of a flexible polymer 
and a semiflexible polymer, where the former takes a spherical globule morphology and the 
latter takes a toroidal morphology in their compact states. We have chosen the inverse 
velocity as A = 10 4 r/a in Fig.l. It can be seen that the pathway of structural transition in 
the loading of a semiflexible polymer is clearly different from that in unloading. 

The force response in a flexible polymer increases at z < 20, maintains a constant value 
at 20 < z < 80, and then increases at z > 80. Correspondingly, a globular flexible polymer 
changes to an ellipsoidal shape and then to a phase-segregated state consisting of a coil and 
a globule, followed by a coil state. This behavior is observed in both loading and unloading, 
and thus a flexible polymer does not exhibit hysteresis. The force response observed in 



22 



Fig^3(b) is called a force plateau, and was predicted by Halperin et al.|2ll. 

In the loading of a semiflexible polymer, a toroidal semiflexible polymer shows phase 
segregation of a coil and a toroid part. In this state, the monomer density decreases in a 
stepwise manner, and correspondingly the force response exhibits a stick- and-release pattern, 
which has been observed in experiments with DNA. At z ~ 140, the polymer undergoes 
a transition to a coiled state. In the unloading process, on the other hand, the phase- 
separated state of a rod and a coil appears at z ~ 125, and the polymer becomes a single 
phase consisting of a rod. At the final stage, a rod becomes more folded or makes a transition 
to a toroid. The force response in a rod-coil phase-segregated state exhibits a force plateau 
and is largely different from that in a toroid-coil state. 
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Next, we consider the cycle in which the initial and final states are coiled states. As 
shown in Fig|21 a semiflexible polymer exhibits a large hysteresis loop, in contrast to a 
flexible polymer. To see the difference more quantitatively, we calculate the hysteresis area 
at various operating speeds. Figure E3 shows the hysteresis area of the force response, i.e. the 
difference in dissipative work between loading and unloading processes. In a flexible polymer, 
hysteresis rapidly decreases as the operating speed decreases. In contrast, hysteresis in a 
semiflexible polymer decreases rather slowly. Large hysteresis remains even at A ~ lOV/a, 
although a semiflexible polymer in bulk makes a transition to a folded state on the time 
scale Afoid ~ 2 x 10 2 r/a. Therefore, we conclude that a semiflexible polymer maintains large 
hysteresis even at an extremely slow operation speed. 

The appearance of notable hysteresis in a semiflexible polymer is associated with the 
significant difference in the pathways between the loading and unloading processes. In un- 
loading, a rod-like structure appears frequently, whereas in loading, a toroidal loop structure 
tends to be generated. As shown in FigsHand|2j the force response strongly depends on the 
structure of the folded part. This structural difference is the origin of the large hysteresis in 
a semiflexible polymer. In contrast, there is almost no difference in the transition pathways 
of a flexible polymer between the loading and unloading processes (FigJ2P). 

Obviously, the hysteresis should disappear at the infinitesimally slow operation where the 
transition between a toroid and a rod is achieved. However in a practical sense, it is difficult 
to observe the transition because of a large free energy barrier between these states. 

IV. SCALING ANALYSIS 

In this section, we will discuss physical meanings of the hysteresis appearing in the cycle. 
Energy balance can be written with the work dW per unit time by the external perturbation 

as, 



where the first and second terms in the r.h.s. are free energy change and dissipative heat, 
respectively. The hysteresis area A in the present system corresponds to the difference be- 
tween total works in loading and unloading process, A = |Wi oa d| — I ^unload I- The dissipative 
heat is equivalent to the entropy production of solvent molecules, 



dW = dF + dQ 



(7) 



dQ = dS. 



(8) 
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In the estimation of the entropy production, only the Stokes drag is considered because the 
hydrodynamic interaction is neglected in the present work. 

We consider a polymer which consists of coiled and collapsed parts. The coiled part has 
n unfolded monomers, and is I in size in the direction parallel to the external force and R in 
the perpendicular direction. The collapsed part is either globule, toroid or rod state. The 
main contribution to the entropy production comes from the coiled part and is 

dS .dR.r, 



where r] is the viscosity. [23J] The monomer size a is hereafter set to be unity. Assuming that 
n monomers behave as a Gaussian statistics in the perpendicular direction to the external 
force, i.e. R ~ n 1//2 , we obtain 

!~" i2 - < io > 

In a flexible polymer, we can evaluate the free energy change to be, 



dF . I 

~ e 

dt K n 



(11) 



The free energy change is reversible between loading and unloading process, and therefore 
it does not contribute to the hysteresis. As a result, the hysteresis in a flexible polymer can 
be expressed as 

A ~ rjzN ~ (12) 

where tr is the Rouse relaxation time of a polymer. 

In a semiflexible polymer, the collapsed parts have different structures between loading 
and unloading process, and therefore the free energy change contributes to the hysteresis. 
The free energy in a loading and unloading process can be written as the summation of three 
contributions; the surface, bending and volume free energy. 

N — n 

^load ~ ^dir 1 + l p 2 e(N-n), (13) 

r i 

d 3 

Unload ~ Hd 2 r 2 + Ip-j^- e(N - n), (14) 

where \i is the surface free energy density. A toroid is regarded as a ring with the radius r\ 
and the width d\, whereas a rod is characterized by the length r 2 and the width d 2 - From 
the conservation of the volume of a collapsed part, the relations dx ~ ((N — n)/ri) l l 2 and 
d 2 ~ ((N — n)/r 2 ) l l 2 have to be satisfied. At a given n, we may assume that the collapsed 
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part rapidly thermalizes, and therefore the surface and bending terms should be balanced. 
As a result, the free energy in a loading and unloading process is given by 

Fioad ~ ll'\N - n) 3 V /5 - <N - n), (15) 
Unload ~ l l J\N - n) 1 ' 2 ^' 2 - e(N - n), (16) 

The entropy production in a semiflexible polymer also arises from the coiled part. Ap- 
proximating that the coiled part consists of cylinders with length l p , we estimate the entropy 
production as 



dS n rjlp ^dR^ 2 
dt l p In l p dt 



(17) 



which is valid l p 3> 1 . The bending energy at the length scale n is of the order of the thermal 
energy k^T so that l p (R/n) 2 ~ 1. As a result, we obtain 

dS rjnz 2 



dt lp In lp ' 



(18) 



From (|15j). ()16|) and (|18j). the hysteresis area in a semiflexible polymer can be estimated as 

A ~ AF+ S> < 19 > 

AF = !j' 5 JV 3 'V 4/5 - ij^iV'V' 2 - (20) 

When £ is small, i.e., the operating speed is small (but still large such that the probability 
of the transition between a rod and a toroid is almost zero), the difference between the 
free energy change in loading and unloading process is larger than the entropy production. 
Therefore, the hysteresis is almost insensitive to the operating speed. The above results 
are qualitatively consistent with those of our simulations (FigOJJ). As shown in (|12p. the 
hysteresis area in a flexible polymer decreases linearly as A = i _1 increases. On the other 
hand, ()19j) indicates that the hysteresis in a semiflexible polymer is almost insensitive to A 
when A is large. 

In this estimation, we neglect the contribution of the relaxation from a rod to a toroid as 
well as from a rod to a more folded rod. However, this process occurs only at small values 
of z and hence does not contribute to the hysteresis when N is large. 



V. TWO- VARIABLE MODEL 

Hysteresis is a general property of a first order phase transition. |25( However, the hys- 
teresis described in the preceding section is different from that in the ordinary case. In a 



semiflexible polymer, the kinetic pathway in the loading process is different from that in 
the unloading process. This is because different structures appear during the transition pro- 
cesses. Therefore, one has to consider two metastable states in the kinetics of the folding and 
unfolding transitions of a semiflexible polymer chain. That is, a transition from a folded 
state to an unfolded state occurs via one metastable state whereas the system traverses 
another metastable state in the reverse process. In order to represent this kind of phase 
transition, one has to introduce at least two variables in the kinetic equations. 

In this section, we consider a model system that provides the above mentioned charac- 
teristic features. We do not intend to develop a quantitative theory specific to a semiflexible 
polymer chain, but rather study the hysteretic property from a general point of view based 
on a simple Ginzburg-Landau type approach. 

We mention a recent theoretical study based on the similar idea. Bartolo et al. have 
investigated dynamic response of adhesion complexes by using the multidimensional energy 
landscape. j3] What they have found is that two alternative trajectories are possible de- 
pending on the loading rate. This is different from our concern that the system relaxes to 
the most stable state via two different metastable states when the process is reversed under 
the same loading rate. 

We start with the free energy in terms of two order parameters X and Y . 

F(X, Y) = -^-X 2 - ^-Y 2 - aXY + ^X 4 + ^Y 4 - h(t)X, (21) 

where fi±, /12 and a are positive coefficients and h is a control parameter depending on 
time. For a suitable set of parameters, the free energy has four local minima, two of which 
correspond to folded and unfolded states. Qualitatively, the variables X and Y correspond, 
respectively, to the density and the orientational order of a polymer with a suitably chosen 
origin of the scale. The coupling between them is incorporated into the third term in the 
free energy ()21|) such that increasing the density increases the orientational order and vice 
versa. 

The molecular dynamics simulations in the preceding section were carried out by changing 
the end-to-end distance at a constant speed. That is, the end-to-end distance was changed 
and the force exerted there was measured. To consider the corresponding situation, we regard 
the external parameter h(t) as the end-to-end distance and assume that there is a one-to-one 
correspondence between the density and the force response. Under these circumstances, we 
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evaluate the hysteresis area in X by changing the external parameter. We fix one parameter 
Hi to be H\ = 100 and examine two situations ^\^> ^2 — 0.1 and \i\ — = 100. When 
H2 = 0.1, there are only two stable minima in the free energy and hence the Y variable is 
irrelevant and corresponds to an ordinary flexible polymer. In contrast, it will be shown 
below that the case with /x 2 = 100 corresponds to a semiflexible polymer. 

The kinetic equations for X and Y are assumed to be given, respectively, by 
dX BF 

— = -L x — + Ci = L x {-X z + viX + aY + h{t)) + £1 (22) 
dY 8F 

-^ = ~L2gy + 6 = L 2 (-Y 3 + /x 2 y + Ql) + £ 2 , (23) 

where Li and L 2 are the Onsager coefficients. The Gaussian random forces £1 and £ 2 satisfy 
the fluctuation-dissipation relation, 

< &(*)&(<') >= 2L i k B T8 ij 8{t - t'). (24) 

We carried out numerical simulations for a set of Langevin equations with L\ — L2 = 1 and 
a = 5.0. The external force is changed linearly in time; 

h{t) = ±vt + h , (25) 

with ho an initial value. 

Figure 4 displays the numerical results of (|2*2*j) and (|23jl . When /i 2 = 0.1, there are 
two stable equilibrium solutions at (X, y) ~ (±12, ±5). With change in the external field 
h(t), transitions occur from one state to the other state. This corresponds to the transition 
between a folded state and an elongated state of a flexible polymer. In fact, the hysteresis 
area decreases as the inverse velocity i _1 is increased, as shown in FigJHl When /i 2 = 100, 
two metastable states appear at (X, Y) ~ (±10, =Fl0) apart from the more stable states 
at (X,Y) ~ (±10, ±10). The transition occurs through a different path via one of these 
metastable states depending on an increase or a decrease in the external field. In this case, 
the hysteresis area does not decrease substantially for a slow operation compared to the case 
with yU 2 = 0.1, as shown in FigJ3J This behavior of hysteresis in Figgis indeed consistent 
qualitatively with the results shown in FigOl obtained by molecular dynamics simulations. 

The scale of the horizontal and vertical axes in Figgis much different from that in FigJBl 
This arises from the uncertainty of the relation between the density and the force response 
as mentioned above and from the difficulty of achieving correspondence to the time scale of 
the folding of a polymer in the artificial model system (}2*2*|) and ([23)1 . 

11 



VI. SUMMARY 



In summary, we have studied the force-strain relation of a bead-spring model of a polymer 
and found the large hysteresis in the loading-unloading cycle for a semiflexible polymer. 
What is different from the ordinary hysteresis is that a large hysteresis is observed even 
with a slow operation, which implies that there are different pathways of conformational 
change between folding and unfolding. This is in contrast to a flexible polymer in which the 
hysteresis area disappears almost entirely with a sufficiently slow operation. This means that 
the conformational change of a flexible polymer can be reduced essentially into a one- variable 
problem. 

In the present paper, we have considered two limiting cases of flexible and semiflexible 
polymers. It would be an interesting problem to investigate systematically the change of 
the hysteresis by changing the bending elasticity k. However this is left for a future study. 

Although most biopolymers such as proteins have much more complicated structures, 
they often exhibit the characteristics of a semiflexible polymer. For example, the secondary 
structures such as a-helix and /?-sheet parts are rather stiff and, therefore, permit a locally 
ordered structure. As a result, even if the initial and final states are the same, the transition 
paths would exhibit a large hysteresis loop, as seen in a semiflexible polymer in the present 
paper. Such kinetic properties are expected to play an important role in biomacromolecules. 
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Figures 
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FIG. 1: Local monomer density and force response under strain in flexible (A, B) and semiflexible 
(C, D) polymers. The black and gray lines show loading and unloading processes respectively. 
Typical snapshots of chain conformation are also shown. The scale bar corresponds to 5a. 
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FIG. 2: Force response in flexible (A) and semiflexible (B) polymers. While the transition pathways 
in loading and unloading processes are the same for a flexible polymer (C), for a semiflexible polymer 
the pathway in loading is largely different from that in unloading (D). 
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FIG. 3: Inverse velocity dependence on hysteresis in flexible (closed triangles) and semiflexible 
(open circles) polymers in our simulations. The solid and broken lines are guides to eyes. 
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FIG. 4: Time trajectory of our model calculation in the XY plane for \ii = 0.1 (A) and \i<i = 100 
(B), and schematic representations of the correspondence to our simulation results in Fig. 0(C) 
and (D). 
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FIG. 5: Inverse velocity dependence on hysteresis at \i2 = 0.1 (closed triangles) and 1x2 = 100 
(open circles) obtained from (|22|) and (J2.3j) . 
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